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1.  Introduction 


Numerical  modeling  of  transport  and  diffusion,  or  the  analysis  of  atmospheric  effects  on  the 
performance  of  sensors,  requires  wind  data  with  high  spatial  (10  to  100  meters)  and  temporal 
(2  to  5  minutes)  resolution  across  a  reasonably  large  (10  x  10  km)  domain.  Traditionally,  there 
have  been  two  approaches  to  generate  such  high  resolution  wind  fields,  and  both  rely  on 
computer  numerical  modeling  of  the  atmosphere.  The  engineering  and  atmospheric  sciences 
communities  usually  simulate  the  wind  field  in  an  urban  area  by  extending  their  computational 
fluid  dynamic  (CFD)  or  large  eddy  simulation  (LES)  models  for  laboratory  control 
environmental  conditions.  Those  models  have  advantages  in  that  their  detailed  treatment  of 
complex  surface  morphologies  allow  for  improved  simulations  of  the  flow  field  structure. 
However,  LES  models  are  generally  not  designed  to  accurately  account  for  the  deep  atmospheric 
structures  that  are  the  major  driving  force  for  the  small  scale  flows.  The  computational  expense 
for  this  type  of  model  is  too  large  to  cover  a  reasonably  sized  urban  area.  On  the  other  hand, 
atmospheric  numerical  weather  prediction  (NWP)  models  usually  attempt  to  represent  a  variety 
of  natural  flow  scales  and  phenomena  including  turbulence,  gravity  waves,  deep  convection, 
fronts,  and  long  waves.  Either  by  increasing  the  NWP  model  resolution  directly  or  alternately  by 
nesting  a  large  eddy  simulation  model  within  a  mesoscale  NWP  model,  one  can  provide  mean 
values  and  associated  turbulence  statistics,  which  are  potential  solutions  for  the  requirement. 
Unfortunately,  it  remains  impractical  for  real-time  applications  in  the  foreseeable  future  because 
of  the  limitations  of  available  computer  power.  Besides  the  huge  requirement  of  computational 
resources,  the  predictability  of  the  forecast  model  in  microscale  remains  debatable  and  is  in  need 
of  further  research.  Eor  applications  requiring  a  combination  of  rapid  and  accurate  results,  a 
well-designed  diagnostic  model,  which  is  exercised  at  specific  times  over  limited  domains,  may 
provide  a  solution. 

One  particular  type  of  diagnostic  model  that  satisfies  some  physical  or  dynamical  constraints  is 
explored  in  this  study.  The  theoretical  basis  for  this  type  of  model  was  developed  by  Sasaki 
(1958;  1970)  using  variational  analysis.  The  general  variational  analysis  defines  an  integral 
function  whose  solution  minimizes  the  variance  of  the  difference  between  the  observed  and 
analyzed  variable  values,  subject  to  physical  constraints.  The  physical  constraints  may  be  any 
combination  of  the  mass,  momentum,  and  energy  conservation  equations.  If  the  physical 
constraints  are  all  three  conservation  equations,  the  problem  becomes  effectively  a  three- 
dimensional  variational  data  assimilation  problem.  Obviously,  this  approach  becomes 
increasingly  more  computationally  intensive  as  constraints  are  added  because  one  not  only  has  to 
solve  the  constraint  equations  but  also  has  to  minimize  the  difference  between  the  observations 
and  the  results  from  the  prediction  model.  If,  however,  we  take  mass  conservation  as  the  only 
constraint,  the  problem  is  greatly  simplified.  However,  this  simplification  will  also  cause  the 
model  to  be  unable  to  simulate  the  turbulence  and  surface  layer  in  the  flow  system.  This 
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question  will  be  addressed  in  detail  in  our  parameterization  of  these  features  in  the  future 
research. 

There  have  been  many  developments  of  this  type  of  model  during  the  last  several  decades.  Ratio 
(1996)  has  provided  a  recent  literature  review  of  mass  consistent  modeling.  In  it,  he  points  out 
that  a  mass  consistent  model  can  provide  a  performance  similar  to  that  of  a  full  dynamical 
simulation  model  in  some  specific  applications.  Sherman  (1978)  and  Dickerson  (1978)  have 
applied  mass  conservation  to  an  analysis  of  atmospheric  flows  over  complex  terrain.  Their 
model  uses  a  local  Cartesian  coordinate,  and  the  earth’s  surface  is  approximated  with  a  constant 
step  in  the  vertical  direction.  The  partition  of  flow  into  both  the  vertical  and  horizontal 
directions  is  based  on  a  simple  empirical  coefficient.  Development  of  other  such  models  by 
Davis,  Bunker,  and  Mutschlecner  (1984)  and  Ross,  Smith,  Manins,  and  Fox  (1988)  employs  a 
terrain- following  sigma  z  coordinate  system.  Connell  (1988)  has  done  extensive  testing  on 
Ross’s  model  using  observational  data,  and  her  results  conclude  that  the  model  performed 
adequately  for  the  flow  over  a  single  mountain.  The  above  authors  attempted  to  formulate  the 
flow  partition  coefficient  in  terms  of  the  Froude  number.  Kitada,  Igarashi,  and  Owada  (1986) 
adopted  a  similar  coordinate  system  in  their  model  and  applied  it  to  a  land-sea  breeze  type  wind 
circulation  pattern.  Moussiopoulos  and  Flassak  (1986)  presented  a  faster  numerical  algorithm, 
using  a  vector  computer,  to  compute  the  wind  field  but  did  not  distinguish  between  the  relative 
adjustment  in  the  horizontal  and  vertical  directions.  These  models  have  typically  covered  100  x 
100  km  domains  at  horizontal  resolutions  of  2  to  4  kilometers.  No  approach  to  date  has  included 
the  effect  of  the  detailed  surface  features  (morphology),  such  as  buildings  and  plant  canopies, 
which  become  important  as  the  resolution  is  increased.  Cionco  (1985)  has  demonstrated,  for  two 
dimensional  simulations,  the  effect  and  impact  of  surface  morphology  on  the  behavior  of  the 
near-surface  flow  as  resolution  is  increased. 

The  objective  in  this  project  is  to  develop  a  robust,  high-resolution  three-dimensional  mass 
consistent  diagnostic  wind  model  (3DSTAT)  which  can  be  applied  in  real-time  simulations 
(Mercurio  et  ah,  2001).  This  model  has  to  be  both  reasonably  accurate  in  representing  the  flow 
physics  and  computationally  efficient.  The  basic  framework  of  the  model  is  discussed  in  this 
document.  The  planned  treatment  of  the  surface  roughness  layer  is  described  in  section  2,  the 
numerical  implementation  of  the  model  is  discussed  in  section  3,  and  several  idealized  tests  of 
the  prototype  of  the  model  are  presented  in  section  4.  The  detailed  treatment  of  the  surface  layer, 
including  the  parameterizations  of  the  turbulent  surface  layer,  the  flow  in  plant  canopies,  the 
turbulent  wakes  by  the  bluff  building  blocks,  and  the  atmospheric  stability  effect,  will  be 
implemented  and  described  in  the  future  as  on-going  processes  of  the  model  development. 
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2.  Model  Description 


2.1  Variational  Formulation 

The  model  is  based  on  the  mass  eonservation  prineiple,  whieh  eliminates  the  divergenee  in  a 
flow  fleld.  That  is,  given  a  limited  number  of  observations  or  a  eoarsely  modeled  wind  field 
over  eomplex  terrain,  the  wind  field  is  physieally  interpolated  in  sueh  way  that  mass 
eonservation  is  satisfied.  Mathematieally,  the  problem  is  to  minimize  the  funetional 

E{u,v,w,A)=l  {u-u^Y  +  {v-v^Y'^Pi  (-  +  —  +  -—)  dhJyt/z  (1) 

dx  dy  dz 


in  whieh  x,  y  are  the  horizontal  eoordinates,  z  the  vertieal  eoordinate,  the  initial 

observed  veloeity  eomponents,  u,  v,  w  the  eorreeted  veloeity  eomponents,  A  the  Lagrange 
multiplier,  and  Pi,  Pi  Gauss  preeision  moduli,  whieh  are  the  wind  veetor  partitioning  faetors  in 
the  horizontal  and  vertical  directions,  respectively.  The  Euler-Lagrange  equations  corresponding 
to  equation  (1)  are 


M  =  M  +■ 


V  =  V  +  ■ 


w  =  w  +  - 


1  dA 

ip^  dx 

1  dA 

dy 

1  dA 

■  2PY  dz 


subject  to  the  boundary  conditions 


du  dv  dw 
—  +  —  +  —  =  0 
dx  dy  dz 


A(u  -u^)  =  0,  yl(v-v°)  =  0,  and  A(w-w^)  =  0 


This  corresponds  to  either  setting  A  =  0  (“flow  through”  free  boundaries)  or  requiring  the  normal 
component  of  the  flow  at  the  boundary  to  remain  unchanged  after  the  adjustment. 

Equations  2  can  be  cast  into  an  equation  for  the  Eagrange  multiplier.  A,  in  terms  of  the  initial 
conditions,  by  differentiating  the  equations  for  u,  v,  and  w,  and  substituting  the  results  into  the 
continuity  equation  to  give  a  Poisson  equation  4.  The  Pi  and  Pi  values  are  assumed  to  be 
constants  throughout  the  small  domain. 
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(4) 


d^A  ^ 

dx^  dy^ 


■  + 


A 

vAy 


d^A 

dz^ 


=  -2  A 


^du°  dv^  dw°^ 
+ - + 


dx  dy  dz 


Without  altering  their  physical  meaning,  let  a  =  (P\/fh)  and  /?i  =  1  so  that  cr  represents  the 
adjustment  of  the  vertical  component  relative  to  the  horizontal  components  (equation  5). 
Determination  of  the  crvalue  is  a  complex  operation  and  it  is  discussed  in  detail  in  section  2.4. 


au 

dx^ 


+ 


d^A  2 
dy^  dz^ 


du^ 
,  dx 


+ 


dw 


o\ 


dz 


(5) 


The  A  value  in  equation  5  can  be  solved  numerically  by  setting  the  boundary  conditions  on  all 
facets  of  the  computation  domain.  The  u,  v,  w  wind  components  then  can  be  computed  from 
equation  2  with  the  A  value  solved  from  equation  5.  The  iterative  convergence  will  be  the  high 
resolution  diagnostic  solution  for  u,  v,  w  for  the  given  boundary  and  coarse  initial  conditions 
(observations).  At  the  lateral  boundaries,  A  is  set  equal  to  zero  to  allow  “flow  through”  in  the 
flow  adjustment.  At  the  bottom  of  the  domain,  “no-flow-through”  conditions  can  be  satisfied  by 
having  the  normal  derivative  vanish,  i.e.,  <9A/<3z  =  0. 

2,2  Vertical  Coordinate 

In  this  model,  a  local  Cartesian  coordinate  {x,y,z)  has  been  chosen  for  the  development  and 
testing.  The  reason  is  that  the  objective  of  the  model  is  to  simulate  microscale  wind  flows  (10  to 
100  m)  in  complex  terrain  or  around  buildings  where  the  computational  domain  and  flow  scale 
are  small.  In  the  terrain-following  coordinate  system,  the  transformation  Jacobian  matrix 
requires  the  derivative  between  the  coordinate  variables  to  exist  in  a  mathematically  continuous 
(smooth)  surface.  The  discontinuity  of  sharp  corners  or  vertical  walls  will  cause  a  singularity  in 
some  element  of  the  Jacobian  matrix.  The  sigma  z  terrain-following  coordinate  must  limit  the 
slope  to  be  less  than  one  vertical  grid  increment  per  horizontal  increment  and  cannot  represent 
severe  terrain  or  city  buildings.  In  this  prototype  of  the  model,  we  use  a  simple  Cartesian 
coordinate.  A  variable  step  vertical  coordinate  following  Tripoli,  Williamson,  and  Garvey 
(2003)  will  be  implemented  in  the  future  to  better  represent  the  topography.  This  vertical 
coordinate  system  represents  a  lower  boundary  by  setting  the  first  layer  depth  to  be  variable  but 
regularly  spaced  above.  The  terrain  surface  will  be  described  almost  exactly  as  prescribed  by  the 
model  horizontal  resolution,  limited  by  the  resolution  of  the  input  terrain  data.  The  idea  of  the 
variable  vertical  step  grid  is  illustrated  in  figure  1.  Consider  the  top  of  the  first  building  block  1, 
where  the  top  of  the  building  is  a  little  greater  than  Z  =  2  but  less  than  Z  =  2.5.  The  constant 
vertical  step  grid  will  approximate  the  building  top  at  Z  =  2,  but  the  variable  step  coordinate  will 
put  the  building  height  at  a  level  indicated  by  placing  a  variable  AZ3  at  this  first  layer. 
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Figure  1.  A  schematic  illustration  of  variable  step  vertical  coordinate. 

2,3  Treatment  of  the  Turbulent  Surface  Layer 

The  wind  near  the  earth’s  surface  is  highly  turbulent  and  chaotic.  Obviously,  the  mass  consistent 
model  will  not  be  able  to  resolve  this  flow.  The  mean  surface  layer  wind  characteristics  have  to 
be  parameterized  according  to  the  observational  data.  For  the  flow  in  and  just  above  plant 
canopies,  the  vertical  wind  profde  proposed  by  Cionco  (1965;  1985)  will  be  used.  The  mean 
wind  profde  model  has  been  verified  through  various  field  experiments  (Cionco,  1985;  1999). 
The  simplicity  of  this  relationship  allows  a  straightforward  coupling  to  the  3DSTAT  model.  The 
turbulence  parameters,  such  as  the  momentum  and  heat  fluxes  and  the  standard  deviation  of  wind 
speed  variation,  need  to  be  parameterized  according  to  the  various  field  measurements.  The  flow 
around  and  above  the  building  blocks  in  an  urban  area  also  requires  special  treatment.  There  are 
some  wind  tunnel  data  (Bossert,  Linn,  Reisner,  Smith,  and  Winterkamp,  1998;  Kastner-Klein, 
Rotach,  Borwn,  Fedorovich,  and  Lawson,  2000;  Brown,  Lawson,  Descroiz,  and  Lee,  2000) 
showing  the  turbulent  wakes  and  boundary  separations  at  the  roof  top  and  sidewalls  of  the 
buildings.  Scarce  data  are  available  for  actual  atmospheric  conditions.  Planned  experiments 
such  as  Joint  Urban  (2003)  will  supply  some  data  for  realistic  modeling  of  the  urban  wind  flow. 
Since  coupling  of  the  turbulent  surface  layer  flow  with  the  upper  level  flow  is  a  complex  topic,  a 
lot  of  detailed  research  and  study  have  to  be  done  at  this  stage.  The  idea  of  parameterizing  the 
flow  in  regions  where  the  model  solution  is  difficult  to  achieve  (or  requires  very  fine  resolution) 
is  not  novel.  This  treatment  is  found  in  many  models,  such  as  boundary  layer  parameterizations 
in  numerical  weather  prediction  models  or  logarithmic  profile  parameterizations  in  some  CFD 
simulations. 
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2.4  Partitioning  the  Flow  Using  the  Richardson  Number 

The  a  value  in  equation  5  defines  the  fiow  partitioning  between  horizontal  and  vertical  motions. 
The  value  of  a  governs  the  adjustment  made  in  the  vertical  component  relative  to  the  horizontal 
components.  In  early  mass-consistent  models,  the  value  was  assumed  to  be  a  small  constant  over 
the  entire  domain  (Sherman,  1978;  Davis  et  al.,  1984)  and  the  value  was  determined  subjectively. 
Snyder,  Thompson,  Eskridge,  Lawson,  Castro,  Lee,  Hunt,  and  Ogawa  (1985)  and  Ross,  Smith, 
Manins,  and  Lox  (1988)  have  related  the  a  value  with  a  hill  Lroude  number  to  account  for  the 
atmospheric  stability.  Recently,  Mercurio,  Williamson,  Chang,  Huynh,  Cionco,  and  Garvey 
(2001)  have  proposed  relating  the  fiow  separation  to  the  bulk  Richardson  number.  This 
partitioning  argument  probably  makes  more  sense,  considering  that  the  model  proposed  here  is 
designed  to  simulate  the  wind  at  an  even  smaller  scale  than  that  of  the  single  hill  treated  by 
Snyder,  Thompson,  Eskridge,  Lawson,  Castro,  Lee,  Hunt,  and  Ogawa,  (1985)  and  Ross,  Smith, 
Manins,  and  Lox  (1988).  Also,  the  a  value  will  be  larger  as  the  flow  scale  becomes  smaller, 
since  the  vertical  motion  in  the  microscale  fiow  can  become  comparable  to  the  horizontal 
motion.  We  defer  this  subject  to  the  next  stage  of  the  model  development. 


3.  Numerical  Methods 


3.1  Finite  Difference  Discretization 


Equation  5  is  discretized  by  a  standard  7-point  method  with  second  order  accuracy  in  all  three 
spatial  directions: 


(Ax)^ 


+  ■ 


(Ay)^ 


-i-ct 


(Az)^ 


=  -2/o  ,  (6) 


in  which is  the  divergence  of  the  wind  field  in  the  last  iteration  (or  the  initial  divergence  in 
case  of  the  first  iteration),  using  a  central  difference  in  space 
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The  “flow-through”  boundary  conditions  are  set  by  letting  A  =  0,  while  the  “no-flow-through” 
boundary  conditions  at  the  terrain  surface  and  building  walls  set  the  normal  derivative  at  the 
point  to  zero  with  a  three-point  forward  (equation  8)  and  backward  difference  (equation  9)  with 
second  order  accuracy. 
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Depending  on  the  direetion  of  the  normal  to  the  boundary,  the  appropriate  normal  direetion  is 
denoted  by  n,  and  the  corresponding  index  is  m.  The  adjusted  wind  field  is  then  computed  from 
the  discretized  equations  of  (2)  with  the  a  value: 
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The  convergence  to  the  final  values  of  u,  v,  w  is  achieved  by  the  iteration  of  equations  6  to  10 
with  a  prescribed  error  tolerance.  The  largest  computational  load  (99%)  arises  in  solving  the 
three-dimensional  Poisson  equation  to  determine  the  X  values.  The  multigrid  method  described 
in  the  next  sub-section  is  employed  to  accelerate  the  computation  of  the  Poisson  equation. 

3,2  Multigrid  Method 

The  discretization  of  the  Poisson  equation  in  a  Cartesian  grid  leads  to  a  system  of  linear 
difference  equations.  This  equation  set  is  non-symmetric,  diagonally  dominant,  and  locally 
dependent  on  the  terrain.  Since  the  equation  set  is  a  very  large  system  (e.g.,  129  x  129  matrix 
for  a  129  X  129  x  129  grid),  an  iterative  method  must  be  applied  to  solve  the  equation  set.  The 
over-relaxation  iteration  method  was  initially  applied  for  this  purpose.  However,  the  over¬ 
relaxation  method  is  quite  slow  to  converge.  Multigrid  algorithms  are  effective  and  fast  in  the 
solution  of  elliptic  equations,  and  they  have  recently  found  many  other  applications.  In  many 
cases,  the  multigrid  method  is  considered  as  the  optimal  method  to  solve  the  elliptical  differential 
equation. 

In  the  tutorial  by  Briggs,  Henson,  and  McCormick  (2000),  the  multigrid  method  was  introduced 
as  a  means  to  accelerate  the  convergence  of  the  relaxation  procedure.  Many  relaxation  schemes 
have  a  smoothing  property,  so  that  the  high  frequency  oscillatory  modes  of  the  error  are 
eliminated  effectively,  but  lower  frequency  smooth  modes  are  dampened  only  slowly.  The 
smoothness  of  the  error  is  relative  to  the  computational  grid  size,  and  a  smooth  mode  in  a  fine 
grid  appears  to  be  high  frequency  to  a  coarse  grid.  The  multigrid  method  takes  advantage  of  this 
property  to  accelerate  the  convergence  by  dealing  with  the  low  frequency  error  in  coarser  grids. 
The  multigrid  method  uses  coarser  grids  recursively  to  relax  the  smooth  mode  error  and 
interpolate  back  to  the  finer  grid.  Coarse  grids  are  also  used  to  compute  an  improved  initial 
guess  for  the  fine  grid  relaxation.  The  full  multigrid  method  (FMG)  applied  in  this  numerical 
model  can  be  described  as  follows: 
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In  order  to  describe  the  FMG  algorithm  clearly,  a  two-grid  correction  scheme  is  described,  first 
as  one-dimensional.  The  three-dimensional  multigrid  method  is  then  a  simple  extension  of  the 
one  dimensional  problem.  The  FMG  V-cycle  scheme  is  used  here  to  solve  the  finite  difference 
equation  LA.  =/at  grid  size  h=  Ax,  where  T  is  a  one-dimension  Poisson  operator  (second 
derivative).  An  appropriate  two-grid  correction  scheme  is 

1 .  Relax  Yi  times  on  l!'  A!'  =f^  on  fine  grid  {Ax  =  h)  with  a  given  initial  guess  f/. 

2.  Compute  the  fine  grid  residual  /  =  f''-Ll'ju’'  and  average  it  to  the  coarse  grid. 

3.  Solve  the  residual  equation  L  e  =  r  for  e  on  coarse  grid  (Ax  =  2h). 

4.  Interpolate  the  coarse  grid  error  to  the  fine  grid  and  then  correct  the  fine  grid 
approximation.  //*<—//*+  e* . 

5.  Relax  72  times  on  =/*  on  the  fine  grid  with  initial  guess 

With  the  two-grid  method  in  mind,  it  is  a  short  step  to  the  multigrid  method.  Instead  of  solving 
the  coarse  grid  residual  equation  exactly,  we  can  get  an  approximate  solution  of  it  by  introducing 
an  even  coarser  grid  and  using  the  two-grid  method.  This  idea  can  be  applied  recursively  down 
to  some  coarsest  grid,  where  the  solution  of  the  error  can  be  found  easily  by  direct  matrix 
inversion  or  iteration.  The  iteration  of  a  multigrid  method  from  finest  grid  to  coarsest  grid  and 
back  to  finest  grid  again  is  called  a  cycle.  If  the  two-grid  iteration  at  each  intermediate  grid  is 
executed  once  only,  it  is  called  a  V  cycle. 

The  FMG  V-cycle  is  a  further  extension  of  the  V  cycle  described  before.  Instead  of  starting  with 
an  arbitrary  approximation  on  the  finest  grid,  which  is  usually  a  poor  guess,  we  start  to  solve  the 
liner  system  at  the  coarsest  grid  and  use  that  solution  to  provide  a  better  guess  for  the  initial  field 
on  the  next  finer  grid.  This  operation  is  conducted  to  the  finest  grid  for  a  good  starting  guess  of 
the  solution  at  each  grid  level.  This  nested  iteration  to  find  the  initial  guess  can  greatly  increase 
the  efficiency  of  the  multigrid  method. 

The  interpolation  method  applied  in  the  three-dimensional  FMG  method  is  a  simple  trilinear 
interpolation.  The  corresponding  averaging  method  is  a  fully  weighted  average.  With  three- 
dimensional  stencil  notation,  the  trilinear  interpolation  (operator  P)  from  coarse  grid  to  fine  grid 
can  be  described  in  the  following  equations  (Wesseling,  1992): 
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in  which  e\,  e2,  and  e3  are  the  unit  directional  vectors,  defined  as 

el  =  (1,0,0),  e2  =  (0,1,0),  e3  =  (0,0,1) 


(12) 


The  restriction  or  averaging  operator  (R)  is  an  adjoint  of  the  operator  P,  whieh  is  a  fully  weighted 
average  in  the  neighboring  points.  The  weighting  factors  for  the  27  points  in  the  three- 
dimensional  steneil  are  expressed  as  three  sliees  of  matrix: 
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in  whieh  is  the  eenter  sliee  that  passes  through  the  eenter  point,  and  i?'*'  and  are  the  first 
and  third  sliees. 
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?(-i) 


The  trilinear  interpolation  and  the  fully  weighted  average  preserve  symmetry  but  are  more 
expensive  eompared  to  other  operators  sueh  as  simple  interpolation  and  simple  injeeting.  A 
detailed  deseription  of  the  interpolation  and  averaging  operators  is  presented  in  Wesseling 
(1992). 

The  relaxation  method  to  solve  the  linear  system  equation  is  the  red-black  Gauss-Seidel  method. 
The  Gauss-Seidel  seheme  is 

=  -T(yv,-/,)  i=i,...,N  (14) 

Ai  7=1 
j*l 


in  whieh  the  N  is  the  number  of  grid  points  in  a  speeifie  grid.  The  grid  points  are  ordered  in  a 
eheeker  board  fashion  in  the  alternating  red  and  blaek  order,  and  the  revising  of  the  red  point 
only  uses  values  at  the  blaek  point,  and  viee  versa.  The  red  and  blaek  points  are  eompletely 
deeoupled  and  the  revising  of  both  groups  may  be  done  in  parallel.  This  property  is  eertainly 
benefieial  on  parallel  eomputers  (Press,  Flannery,  Teukolsky,  and  Vetterling,  1992).  The 
eonvergenee  theory  for  the  multigrid  method  is  beyond  the  seope  of  this  report,  and  the  reader  is 
referred  to  Wesseling  (1992)  and  Briggs,  Henson,  and  MeCormiek  (2000). 

The  interpolation,  averaging,  and  relaxation  operator  use  exaetly  the  same  subprogram  at  every 
grid  level.  The  numerieal  routines  are  implemented  via  a  standard  FORTRAN  (Formula 
Translator)  90  eode  strueture.  The  dynamieal  memory  alloeation  and  the  reeursive  subroutine 
are  used  to  improve  the  effieieney  of  the  numerieal  eode. 
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4.  Initial  Test  Results  of  the  Model 


4,1  Test  With  the  Potential  Flow  Solutions 

The  model  formulation  and  eoding  require  very  rigorous  testing  to  detect  potential  coding  errors. 
A  logical  first  test  is  to  compare  with  analytical  solutions  if  available.  The  model  formulation  of 
the  minimization  of  the  functions  resulted  in  a  Poisson  equation  5  with  respect  to  the  Lagrange 
multiplier  X.  By  letting  a  =  1  (the  wind  is  adjusted  equally  in  the  vertical  and  horizontal 
directions),  equation  5  can  rewritten  in  vector  form; 

V^A  =  -2V*fo  (15) 

The  X  represents  a  velocity  potential  for  non-viscous,  irrotational  flow  if  the  initial  wind 
divergence  V  •Vg  =0.  This  condition  is  satisfied  for  the  uniform  background  wind.  Thus,  a 

mass  consistent  numerical  model  simulation  of  an  initially  uniform  wind  flow  over  a  terrain  is 
equivalent  to  the  potential  flow  solution.  The  potential  flow  can  be  expressed  in  a  more  general 
sense  (Ross,  Smith,  Manins,  and  Fox,  1988):  If  the  initial  wind  Vg  =  (uo,  vo,  wo)  can  be  separated 
as 


Uo  =  uo{x),  Vo  =  vo(y)  and  ivo  =  wo{z) 


(16) 


then  the  final  wind  satisfies  a  velocity  potential  O 


do 


(17) 


in  which 


o  = 


J  Ugdx  +  I  Vgdy  +  I  Wgdz 


A 

+  — 
2 


(18) 


We  then  have  an  equation  that  is  identical  with  (15) 


V"O  =  V»F„+-VU  =  0 
°  2 


(19) 


There  are  some  analytical  solutions  for  the  potential  flow  around  simple  two-dimensional  or 
three-dimensional  geometrical  objects  such  as  a  cylinder,  sphere,  or  ellipse  (Milne-Thompson, 
1960).  However,  for  complex  geometrical  obstacles,  analytical  solutions  are  not  possible  and 
numerical  simulation  is  the  preferred  method.  Nevertheless,  the  simple  geometry  analytical 
solutions  are  useful  for  the  testing  of  the  numerical  model.  The  following  are  some  model  tests 
compared  to  some  simple  potential  flow  analytical  solutions.  The  potential  flow  theory  is  a  good 
approximation  for  the  flow  over  an  obstacle  in  the  neutral  condition  and  when  the  Froude 
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number  {Fr)  is  much  greater  than  1 ,  except  in  the  region  of  turbulent  wake  in  the  lee  side  of  the 
flow.  The  Froude  number  is  defined  as 


F 


Nh 


(20) 


in  which  is  the  uniform  wind  speed  approaching  the  hill,  h  is  the  height  of  the  hill,  and 

N  =  -  —  ,  the  Brunt- Vaisala  frequency.  The  condition  of  »  1  can  be  found  in  the 

y  p  dz 

condition  when  (1)  there  is  a  homogenous  and  well-mixed  condition  (F  ^  0)  in  which  the 


Reynolds  number  is  usually  large  (Baines,  1995),  or  (2)  a  small  hill  height  (h),  given  the  flow  is 
stratified.  The  parameterization  of  the  lee  side  turbulence  will  be  studied  in  detail  in  the  next 
stage  of  the  research.  The  potential  flow  theory  is  to  be  compared  here  only  for  the  purpose  of 
checking  the  computational  methods  and  the  flow  upstream  and  above  the  turbulence  wake. 


a.  Potential  Flow  Over  a  Flemisphere 


A  simple  hemisphere  with  a  radius  of  2.5  km  is  used  for  this  case.  The  potential  flow  analytical 
solution  for  this  case  is  an  axis  symmetric  flow  in  a  polar  coordinate  {r,6)  (Milne-Thompson, 
1960); 


Fr  =  F„  cos  6{l  -  -^) 
a 

V,  =  -V^  sin  d(l  +  ^) 

2a 


(21) 


in  which  and  Vg  are  the  wind  velocity  components  normal  and  tangential  to  the  hemisphere 
surface,  respectively,  9  is  the  elevation  angle  of  the  line  normal  to  the  point,  is  the  uniform 
prescribed  wind  speed,  a  is  the  distance  from  the  hemisphere  center,  and  r  is  the  radius  of  the 
hemisphere.  The  potential  solution  is  simulated  with  the  3DSTAT  model  on  a  10-km  x  10-km  x 
10-km  domain.  The  grid  number  is  129  x  129  x  129  and  dx  =  dy  =  dz  =  77.5m.  The  radius  of 
the  hemisphere  is  2.5  km  and  the  uniform  wind  speed  is  10  m/s.  The  comparison  of  the 
analytical  solution  and  the  numerical  solution  is  plotted  in  figure  2.  The  left  panel  is  the 
analytical  solution  and  the  right  panel  is  the  3DSTAT  simulation.  In  both  solutions,  the  top 
panel  displays  the  horizontal  cross  section  at  z  =  1.55  km  and  the  bottom  panels  shows  the 
vertical  cross  section  at  y  =  5  km.  Overall,  the  numerical  simulation  captured  the  potential  flow 
solution  and  showed  the  wind  turning,  climbing,  accelerating  on  the  mountain  top,  and 
approaching  zero  at  the  stagnation  point. 
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Analytical  solution 
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Figure  2.  Comparison  of  the  potential  flow  analytieal  solution  (left  panel  figures)  with  the 
numerieal  model  solution  (right  panel  figures).  (The  top  panels  are  the  horizontal 
eross  sections  at  z  =  1.55  km.  The  bottom  panels  are  the  vertical  cross  sections  at 
y  =  5  km.) 

b.  Potential  Flow  Over  Two-Dimensional  Ridges 

There  are  potential  flow  solutions  for  many  two-dimensional  geometries.  Two  have  been  chosen 
for  the  purpose  of  testing  the  model.  The  flow  over  the  half  cylinder  and  the  flow  over  a  two- 
dimensional  “Witch-of- Agnesi”  (representing  long  ridges  in  the  y  direction),  are  presented  here. 
An  important  feature  for  stratified  flow  over  those  types  of  ridges  (Queney,  1948;  Long,  1953) 
and  over  a  three-dimensional  Witch-of- Agnesi  (Smith,  1989)  is  that  solutions  based  on  linear 
perturbation  theory  show  gravity  waves  generating  along  the  lee  side  of  the  ridges  for  large 
Froude  numbers.  The  diagnostic  model  presented  here  will  not  show  this  type  of  waves  because 
it  does  not  address  momentum  conservation.  Flowever,  for  a  small  obstacle  or  for  the  well- 
mixed  atmospheric  condition,  the  potential  flow  is  a  valid  approximation  except  in  the  a  flow 
separation  region  such  as  lee  wake.  Figure  3  compares  the  potential  flow  analytical  solutions 
(Milne-Thompson,  1960;  Queney,  1948)  and  the  3DSTAT  simulation.  The  top  panel  is  for  the 
flow  over  a  half  cylinder  and  the  bottom  panel  for  the  flow  over  the  two-dimensional  Witch-of- 
Agnesi  ridge.  The  streamlines  in  both  simulations  match  the  analytical  solution  reasonably  well. 
However,  the  model  simulations  have  some  differences  near  the  surface  of  these  ridges. 
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Figure  3.  Streamlines  of  flow  around  infinite  long  half  cylinder  ridge  (top  panel);  streamlines 
of  flow  over  a  2-D  Witch-of-Agnesi  ridge  (bottom  panel). 

4,2  Flow  Around  and  Over  Building  Blocks 

Wind  flow  around  building  blocks  or  building  canyons  is  a  very  complex  problem.  There  have 
been  a  number  of  wind  tunnel  studies  and  reviews  of  wind  flow  around  simple  obstacles  (Hunt, 
Abell,  Peterka,  and  Woo,  1978)  or  groups  of  building  blocks  (Bossert,  Linn,  Reisner,  Smith,  and 
Winterkamp  1998;  Kastner-Klein,  Rotach,  Brown,  Fedorovich,  and  Lawson  2000;  Brown, 
Lawson,  Descroiz,  and  Lee,  2000).  For  an  isolated  cubical  obstacle.  Hunt,  Abell,  Peterka,  and 


13 


Woo,  (1978)  showed  that  several  areas  in  the  immediate  vicinity  of  the  obstacle  have  re¬ 
circulations  and  higher  turbulent  intensities  because  of  the  steep  sides  and  sharp  corners 
(figure  4).  On  the  central  plane  of  symmetry  of  fluid  flow  past  a  cuboid  there  are  the  upstream 
re-circulation  region,  the  roof-separated  flow  region,  and  the  main  wake  with  reversed  flow  on 
the  lee  side.  At  the  side  walls,  there  are  also  flow  separations.  Brown,  Lawson,  Descroiz,  and 
Lee  (2000)  showed  that  when  a  wind  flowed  perpendicularly  over  a  seven-building  array,  the 
flow  separated  at  the  top  of  the  first  building  but  not  over  the  remaining  downstream  buildings. 
The  flow  at  the  surface  was  nearly  identical  in  all  street  canyons,  where  there  was  a  re¬ 
circulation  region  with  a  strong  reversed  flow  at  0.6  building  height. 


Figure  4.  A  schematic  graph  shows  the  mean  flow  pattern  over  and  around  a  wall-mounted  box. 

Figure  5  shows  the  3DSTAT  simulation  (top  panel)  of  the  flow  over  three  building  blocks 
compared  with  the  wind  tunnel  results.  Obviously,  we  do  not  expect  the  3DSTAT  model  to 
simulate  the  canyon  vortex  and  the  boundary  layer  separation  at  the  top  of  the  first  building 
because  the  model  does  not  contain  the  momentum  conservation.  The  flow  outside  these 
regions,  which  can  be  approximated  by  the  potential  flow  theory  (Baines,  1995),  is  similar  to  the 
wind  tunnel  results.  The  wake  area  and  the  flow  separation  area  have  to  be  parameterized  in 
further  model  development.  Figure  6  shows  another  3DSTAT  solution  for  the  wind  flow  around 
and  above  a  group  of  five  buildings.  The  first  two  buildings  on  the  left  (see  the  top  panel)  are 
13  meters  high,  while  the  others  are  8  meters  high.  The  streamlines  show  that  the  3DSTAT 
simulation  captured  the  expected  flow  around  the  buildings  (top  panel)  and  the  building  blocking 
effects.  The  bottom  panel  shows  the  flow  over  the  vertical  cross  section  at  y  =  30m.  The  wind  at 
the  top  building  also  shows  acceleration.  Flowever,  the  wakes  at  the  lee  side  of  the  buildings  are 
not  simulated  because  momentum  conservation  is  not  imposed. 
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Figure  5.  3DSTAT  simulation  of  wind  flow  around  and  over  3-D  building  blocks.  (The  top  panel 
is  a  horizontal  cross  section  at  z  =  5m,  and  the  bottom  panel  is  a  vertical  cross  section  at 
y  =  30m.) 


15 


Mind-tunnel  dat^  (U*5.  EPH,  frciH  Bk»s&ert  et  al^) 

25U 

1&D 
1DD 
bO 
D 

-15D  0  150  3og  450  SDO 

Trrrrrrr 


Figure  6.  Flow  over  a  3-D  block  array.  (3DSTAT  simulation  (top  panel)  and  the  wind  tunnel  data  (bottom 
panel,  U.S.  EPA)). 

4,3  Computational  Efficiency  of  the  Model 

The  multigrid  method  applied  in  the  computation  very  much  improves  the  computational 
efficiency.  Before  implementation  of  the  multigrid  method,  we  developed  a  test  version  of  the 
model  use  a  simple  red-black  (RB)  Gauss-Seidel  over-relaxation  method  to  solve  the  Poisson 
equation.  Table  1  lists  the  control  processing  unit  (CPU)  times  required  for  the  test  cases  with 
both  versions  of  the  model.  The  testing  was  done  on  a  two-processor  2-GHz  Pentium  4  Linux 
Dell  computer.  Basically,  this  version  of  the  model  takes  about  20  to  30  times  more  CPU  time  to 
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run,  compared  with  the  latest  model  version  implemented  with  the  multigrid  model  (depending 
on  the  complexity  of  the  underlying  terrain).  Another  advantage  of  the  multigrid  method  is  that 
the  eomputation  time  is  not  proportionally  inereased  as  the  terrain  surface  beeomes  more 
eomplex.  Indeed,  test  case  2  over  case  4  only  showed  about  a  1.5-minute  CPU  time  inerease  in 
the  multigrid  version  model,  while  it  had  a  2-hour  inerease  of  CPU  time  with  the  simple  RB 
Gauss-Seidel  over-relaxation.  We  have  tested  this  eode  for  sealability  for  multiproeessor 
eomputers  using  a  simple  OpenMP  compiling  option.  The  results  showed  that  the  eode  scales 
fairly  well  up  to  four  proeessors  with  about  a  1.8  time  aceeleration.  Indeed,  the  RB  Gauss-Seidel 
relaxation  has  very  good  properties  when  the  red  and  blaek  points  are  deeoupled  so  that  parallel 
proeessing  ean  be  easily  applied. 


Table  1.  Comparison  of  CPU  time  (minutes)  with  different  numerical  methods. 


Test  Cases 

Topographies 

Simple  RB 
Gauss-Seidel 
method 

Multigrid  method 

1 

Hemisphere  (3D) 

124 

6 

2 

Half-cylinder  (2D) 

119 

6.4 

3 

Witch-of- Agnesi  (2D) 

120 

7.1 

4 

Building  array  (2D) 

218 

7.5 

5 

Five  buildings  (3D) 

207 

7.8 

5.  Summary  and  Conclusions 


This  doeument  deseribes  a  prototype  high  resolution,  three-dimensional,  computationally 
efficient,  diagnostic  model  for  flow  over  complex  terrain  with  a  mass-consistent  approaeh.  The 
differenees  between  the  eurrent  model  and  similar  approaehes  are  in  the  vertieal  eoordinate,  the 
lower  boundary  conditions,  and  the  ehoice  of  numerical  method.  This  model  will  use  a  step 
vertical  coordinate  whieh  is  better  suited  for  very  steep  topography  and  will  include  the  effeets  of 
not  only  topography  but  also  small  surfaee  features  sueh  as  trees  and  buildings  on  the  overall 
flow.  The  numerical  implementation  takes  advantage  of  a  multigrid  method  which  greatly 
improves  the  eomputation  speed.  The  framework  of  the  model  and  assoeiated  implementations 
have  been  deseribed  here,  several  preliminary  test  eases  for  the  model  have  been  demonstrated, 
and  defieiencies  of  the  model  have  been  addressed.  The  proposed  parameterization  of  the 
turbulent  surface  layer  and  the  flow  partition  will  be  performed  in  our  future  researeh  with  field 
observation  data. 
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